01 February 2018

Why?

Mapping in Ecology?

  1. show where your plots are

  2. show how variables are distributed spatially

  3. spatial analyses

  4. Many, many tools

R as a GIS?

  1. Workflow

  2. Quite efficient

  3. Interface

Workflow

R for everything (almost)

  1. import your data
  2. format your data
  3. analyze your data
  4. visualize your data
  5. export your data
  6. create your own functions/packages

Efficiency

  1. well-defined classes
  2. many formats (+ convert different file)
  3. operations on geometries
  4. statistical analyses
  5. visualization

Interface

Overview

Turning R into a powerful GIS

1- Define classes

  • sp
  • raster
  • sf

2- Topology operations

  • rgeos
  • sf

3- Import / export

  • rgdal
  • sf

Turning R into a powerful GIS

4- Analyses

  • spdep
  • dismo

5- Creating and editing maps

Turning R into a powerful GIS

6- But

  • Geo-referencing
  • Visualizing large spatial objects

Package sp

Classes

Classes / Functions Contents
Points list of points (set of coordinates)
SpatialPoints list of points + CRS
SpatialPointsDataPoints list of points + CRS + attribute table
Line a line (set of coordinates)
Lines list of lines
SpatialLines list of lines + CRS
SpatialLinesDataFrame list of lines + CRS + attribute table

Example: SpatialPointsDataPoints

long lat var1 var2
-80.07779 43.23585 -0.9424900 8.822271
-81.58100 43.29668 0.4211175 9.058452
-81.27616 42.45650 2.1195748 4.783351
-81.73602 43.31621 0.3377636 8.061516
-81.53676 43.28936 0.8181175 4.986830
-80.23726 42.64735 -0.5898903 4.639604

Example: SpatialPointsDataFrame

library(sp)
## Loading required package: methods
mysp <- SpatialPointsDataFrame(
  coords = mydata[,1:2],
  data = mydata[,3:4],
  proj4string = CRS(
    "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0"
  )
)

Example: SpatialPointsDataFrame

class(mysp)
## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"
isS4(mysp)
## [1] TRUE
slotNames(mysp)
## [1] "data"        "coords.nrs"  "coords"      "bbox"        "proj4string"

Example: SpatialPointsDataFrame

mysp@proj4string
## CRS arguments:
##  +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0
head(mysp@data)
##         var1     var2
## 1 -0.9424900 8.822271
## 2  0.4211175 9.058452
## 3  2.1195748 4.783351
## 4  0.3377636 8.061516
## 5  0.8181175 4.986830
## 6 -0.5898903 4.639604

Change projection: spTransform

(mysp2 <- spTransform(mysp, CRS=CRS("+proj=merc +ellps=GRS80")))
##            coordinates        var1      var2
## 1  (-8914219, 5318662) -0.94249004 8.8222714
## 2  (-9081556, 5327927)  0.42111746 9.0584523
## 3  (-9047620, 5200757)  2.11957478 4.7833508
## 4  (-9098812, 5330904)  0.33776361 8.0615156
## 5  (-9076630, 5326812)  0.81811746 4.9868304
## 6  (-8931971, 5229491) -0.58989032 4.6396036
## 7  (-9049053, 5341975) -0.14429222 7.9327105
## 8  (-9082845, 5365187)  0.73672200 4.9249723
## 9  (-9010743, 5223804)  0.24897363 4.0290142
## 10 (-9042984, 5408682) -0.08106143 3.0501115
## 11 (-8992699, 5358424) -0.46859518 0.6339195
## 12 (-8945265, 5158475)  0.70483256 1.5400158
## 13 (-8976821, 5431270) -0.62466071 0.4572709
## 14 (-9095845, 5201313)  1.93421826 0.9320450
## 15 (-9110967, 5255899)  0.77250278 6.2389721
## 16 (-9024592, 5159324) -1.09241480 7.7116822
## 17 (-8935065, 5326269) -0.41269577 4.4730041
## 18 (-9070013, 5264491) -0.32887855 4.7979019
## 19 (-9011305, 5231279) -1.14365387 1.1729072
## 20 (-9039393, 5432599)  1.54523563 7.8279072

Package raster

Classes

0- SpatialGrid / SpatialPixel (package sp)

3- RasterLayer (package raster)

4- RasterStack (package raster)

5- RasterBrick (package raster)

Example:`raster`

library(raster)
## 
## Attaching package: 'raster'
## The following object is masked from 'package:magrittr':
## 
##     extract
ras1 <- raster(matrix(runif(100*100,0,10),ncol=100,nrow=100),
    crs=CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"),
    xmn=-82, xmx=+80, ymn=+42, ymx=+44)

Example:`raster`

ras1
## class       : RasterLayer 
## dimensions  : 100, 100, 10000  (nrow, ncol, ncell)
## resolution  : 1.62, 0.02  (x, y)
## extent      : -82, 80, 42, 44  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : in memory
## names       : layer 
## values      : 0.0001147902, 9.997746  (min, max)

Retrieving free GIS data: getData

head(getData("ISO3"))
##   ISO3        NAME
## 1  ABW       Aruba
## 2  AFG Afghanistan
## 3  AGO      Angola
## 4  AIA    Anguilla
## 5  ALA       Åland
## 6  ALB     Albania

Example:`getData`

## Country level:
mapBEL0 <- getData(name="GADM", country="BEL", path="./assets", level=0)
mapBEL1 <- getData(name="GADM", country="BEL", path="./assets", level=1)
tminW <- getData(name="worldclim", var="tmin", res=10, path="./assets")
mapBEL0
## class       : SpatialPolygonsDataFrame 
## features    : 1 
## extent      : 2.555356, 6.40787, 49.49722, 51.50382  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## variables   : 68
## names       : OBJECTID, ID_0, ISO, NAME_ENGLISH, NAME_ISO, NAME_FAO,              NAME_LOCAL, NAME_OBSOLETE,            NAME_VARIANTS, NAME_NONLATIN, NAME_FRENCH, NAME_SPANISH, NAME_RUSSIAN, NAME_ARABIC, NAME_CHINESE, ... 
## value       :        1,   23, BEL,      Belgium,  BELGIUM,  Belgium, Belgique|Belgie|Belgien,              , Belgique|Belgien|Belgium,              ,   Belgique ,     Bélgica ,     Бельгия ,      بلجيكا,      比利时 , ...

Package rgdal

Drivers

library(rgdal)
## rgdal: version: 1.2-16, (SVN revision 701)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.1.2, released 2016/10/24
##  Path to GDAL shared files: /usr/share/gdal/2.1
##  GDAL binary built with GEOS: TRUE 
##  Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
##  Path to PROJ.4 shared files: (autodetected)
##  Linking to sp version: 1.2-5
  1. writeOGR / writeGDAL: to write spatial objects
  2. readOGR/readGDAL: to read spatial files

Drivers

head(ogrDrivers())
##         name                     long_name write  copy isVector
## 1 AeronavFAA                   Aeronav FAA FALSE FALSE     TRUE
## 2 AmigoCloud                    AmigoCloud  TRUE FALSE     TRUE
## 3     ARCGEN             Arc/Info Generate FALSE FALSE     TRUE
## 4     AVCBin      Arc/Info Binary Coverage FALSE FALSE     TRUE
## 5     AVCE00 Arc/Info E00 (ASCII) Coverage FALSE FALSE     TRUE
## 6        BNA                     Atlas BNA  TRUE FALSE     TRUE

Export

writeOGR(mysp, dsn="./assets", layer="mypoints",
    driver="ESRI Shapefile", overwrite_layer=TRUE)

Import

mysp2 <- readOGR(dsn="assets/", layer="mypoints")
## OGR data source with driver: ESRI Shapefile 
## Source: "assets/", layer: "mypoints"
## with 20 features
## It has 2 fields
## Roads find on http://www.diva-gis.org/Data
canroads <- readOGR(dsn="assets/", layer="CAN_roads")
## OGR data source with driver: ESRI Shapefile 
## Source: "assets/", layer: "CAN_roads"
## with 16915 features
## It has 5 fields

Note: there are also functions to import/export raster in raster

Package rgeos

Load the package

library(rgeos)
## rgeos version: 0.3-26, (SVN revision 560)
##  GEOS runtime version: 3.5.1-CAPI-1.9.1 r4246 
##  Linking to sp version: 1.2-5 
##  Polygon checking: TRUE

Load the package

buf <- gBuffer(mapBEL0, width=0.5)
## Warning in gBuffer(mapBEL0, width = 0.5): Spatial object is not projected;
## GEOS expects planar coordinates
diff <- gDifference(gBuffer(mapBEL0, width=0.5), gBuffer(mapBEL0, width=0.1))
## Warning in gBuffer(mapBEL0, width = 0.5): Spatial object is not projected;
## GEOS expects planar coordinates
## Warning in gBuffer(mapBEL0, width = 0.1): Spatial object is not projected;
## GEOS expects planar coordinates

Package mapview

Import package

library(mapview)

NB: it uses the leaflet package.

Quick examples

mapview(mysp, cex = 'var2')@map

Quick examples

mapview(mapBEL1)@map

More

Editing a map

A very basic map – Shapefile

plot(mapBEL0)

A very basic map – Raster

class(tminW)
## [1] "RasterStack"
## attr(,"package")
## [1] "raster"

A very basic map – Raster

plot(tminW)

Customize a map – Shapefile

plot(mapBEL0, border='grey15', col='#E6E6E6', lwd=1.6)
plot(mapBEL1, lty=2, lwd=0.9, add=T)
points(4.3513, 50.8471, pch=19, col="#27df9d")
text(4.3513, 50.8471, text="Brussel", pos=3)

Customize a map – Shapefile

Customize a map – Shapefile of roads

package sf

Simple features?

Importing spatial object

Doing simple operations

Exporting

Examples

Let’s practice

Installation

Part 2 – Showcases

Data frame with Ontario’s lakes

  • the .csv file below is available in the repository
lakedf <- read.csv('assets/lakeOnt.csv')
head(lakedf)
##   X LAKEID LONGITUDE_DD LATITUDE_DD Area_ha
## 1 1 L32173    -80.03361    49.90417   174.9
## 2 2 L21001    -91.68278    48.68028   130.2
## 3 3 L44001    -84.19667    46.91861   279.6
## 4 4 L25002    -90.43111    48.16639   114.1
## 5 5 L44002    -84.30528    47.06000   137.3
## 6 6 L21003    -91.34528    48.24806  2993.9

Make it a SpatialPointDataFrame

library(sp)
lakesp <- SpatialPointsDataFrame(
  coords = lakedf[,3:4],
  data = lakedf[,c(1,4)],
  proj4string = CRS(
    "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0"
  )
)

Your SpatialPointDataFrame

lakesp
## class       : SpatialPointsDataFrame 
## features    : 546 
## extent      : -95.06111, -76.01194, 43.05722, 54.52056  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0 
## variables   : 2
## names       :   X, LATITUDE_DD 
## min values  :   1,    43.05722 
## max values  : 546,    54.52056

Retrieve altitude for Canada

altCAN <- getData(name="alt", country="CAN", path="./assets/")
altCAN
## class       : RasterLayer 
## dimensions  : 4992, 10620, 53015040  (nrow, ncol, ncell)
## resolution  : 0.008333333, 0.008333333  (x, y)
## extent      : -141.1, -52.6, 41.6, 83.2  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
## data source : /home/kevcaz/Github/Tutorials/mapsWithR/docs/assets/CAN_msk_alt.grd 
## names       : CAN_msk_alt 
## values      : -108, 5800  (min, max)

Retrieve altitude for Canada

Retrieve provinces’ boundaries

bouCAN <- getData(country='CAN', level=1, path="./assets/")
bouCAN
## class       : SpatialPolygonsDataFrame 
## features    : 15 
## extent      : -141.0069, -52.61889, 41.67693, 83.11042  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## variables   : 13
## names       : OBJECTID, ID_0, ISO, NAME_0, ID_1,  NAME_1, HASC_1, CCN_1, CCA_1,     TYPE_1, ENGTYPE_1, NL_NAME_1,                                       VARNAME_1 
## min values  :        1,   42, CAN, Canada,    1, Alberta,  CA.AB,    NA,    10,   Province,  Province,          ,                                                 
## max values  :       15,   42, CAN, Canada,   13,   Yukon,  CA.YT,    NA,    62, Territoire, Territory,          , Yukon Territory|Territoire du Yukon|Yukon|Yuk¢n

Retrieve provinces’ boundaries

bouCAN@data[,1:8]
##    OBJECTID ID_0 ISO NAME_0 ID_1                    NAME_1 HASC_1 CCN_1
## 1         1   42 CAN Canada    1                   Alberta  CA.AB    NA
## 2         2   42 CAN Canada    2          British Columbia  CA.BC    NA
## 3         3   42 CAN Canada    3                  Manitoba  CA.MB    NA
## 4         4   42 CAN Canada    4             New Brunswick  CA.NB    NA
## 5         5   42 CAN Canada    5 Newfoundland and Labrador  CA.NF    NA
## 6         6   42 CAN Canada    6     Northwest Territories  CA.NT    NA
## 7         7   42 CAN Canada    7               Nova Scotia  CA.NS    NA
## 8         8   42 CAN Canada    8                   Nunavut  CA.NU    NA
## 9         9   42 CAN Canada    8                   Nunavut  CA.NU    NA
## 10       10   42 CAN Canada    8                   Nunavut  CA.NU    NA
## 11       11   42 CAN Canada    9                   Ontario  CA.ON    NA
## 12       12   42 CAN Canada   10      Prince Edward Island  CA.PE    NA
## 13       13   42 CAN Canada   11                    Québec  CA.QC    NA
## 14       14   42 CAN Canada   12              Saskatchewan  CA.SK    NA
## 15       15   42 CAN Canada   13                     Yukon  CA.YT    NA

Canada – Ontario

bouCAN[11,]
## class       : SpatialPolygonsDataFrame 
## features    : 1 
## extent      : -95.15598, -74.34605, 41.67693, 56.86972  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## variables   : 13
## names       : OBJECTID, ID_0, ISO, NAME_0, ID_1,  NAME_1, HASC_1, CCN_1, CCA_1,   TYPE_1, ENGTYPE_1, NL_NAME_1,    VARNAME_1 
## value       :       11,   42, CAN, Canada,    9, Ontario,  CA.ON,    NA,    35, Province,  Province,          , Upper Canada

Canada – Ontario

Canada – Ontario

plot(crop(altCAN, bouCAN[11,]@bbox))
par(mar=c(0,0,0,0))

Canada – Ontario

Canada – Ontario elevation

altONT <- rasterize(bouCAN[11,], crop(altCAN, bouCAN[11,]@bbox), mask=TRUE)
mapview(altONT)+mapview(lakesp)

Canada – Ontario elevation

Canada – Ontario elevation

par(mar=c(0,0,0,0))
plot(altONT)
plot(lakesp, add=T)

Canada – Ontario elevation

Resources

Useful links to use R as a GIS

Useful links to get data